load("./data/elk_processed.rda")
library(ggplot2)
ggplot(data=elk_processed, aes(x=lon, y=lat)) +
geom_path(size=0.5, color="darkgrey") +
theme_classic() +
facet_wrap(~id, scale="free", ncol=3)
library(dplyr)
elk_mig <- elk_processed |>
filter(id %in% c("GP2", "YL25", "YL29", "YL73", "YL77", "YL78"))
ggplot(data=elk_mig, aes(x=datetime, y=lat)) +
geom_path(size=0.5, color="darkgrey") +
theme_classic() +
facet_wrap(~id, scale="free", ncol=2)
# Data Regularization (AniMotum)
install from my R-universe repository
need latest TMB and Matrix package versions also
install.packages("aniMotum",
repos = c("https://cloud.r-project.org",
"https://ianjonsen.r-universe.dev"),
dependencies = TRUE)
remotes::install_version('TMB', version = '1.9.10')
install.packages("Matrix")
library(aniMotum)
data should be structured with columns: id, date (POSIXct), lc, lon, lat and (optional) additional locational error information smaj, smin, eor (Argos tags)
lc (location class) can include the following values: 3, 2, 1, 0, A, B, Z, G, or GL. The latter two are for GPS locations and ‘Generic Locations’, respectively. Class Z values are assumed to have the same error variances as class B. By default, class G (GPS) locations are assumed to have error variances 10x smaller than Argos class 3 variances, but unlike Argos error variances the GPS variances are the same for longitude and latitude
elk_mig2 <- elk_mig |>
mutate(lc = "G", date = datetime) |>
dplyr::select(id, date, lc, lon, lat)
check sampling rate
dtime <- function(t, ...) {difftime(t[-1], t[-length(t)], ...) %>% as.numeric}
sample.rate <- elk_mig2 |>
group_by(id) |>
arrange(date) |>
mutate(dtime = c(0,round(dtime(date, units ="hours")))) |>
mutate(mode_dtime = as.numeric(names(sort(table(dtime), decreasing=TRUE))[1])) |>
data.frame() |>
ungroup()
barplot(table(sample.rate$dtime))
use sample rate of 2 hours
vmax: m/s
time.step: hrs
random walk: better for large data gaps than crw
elk_ssm_fit <- fit_ssm(elk_mig2,
vmax = 20,
model = "rw",
time.step = 2,
control = ssm_control(verbose = 0))
summary(elk_ssm_fit)
## Animal id Model Time n.obs n.filt n.fit n.pred n.rr converged AICc
## GP2 rw 2 6082 0 6082 3546 . TRUE 4542.8
## YL25 rw 2 6413 1 6412 4664 . TRUE 14777.5
## YL29 rw 2 5405 3 5402 3806 . TRUE 10070.9
## YL73 rw 2 9614 0 9614 3076 . TRUE 8472.5
## YL77 rw 2 2802 0 2802 3230 . TRUE 13299.3
## YL78 rw 2 3929 0 3929 2810 . TRUE 6841.3
##
## --------------
## GP2
## --------------
## Parameter Estimate Std.Err
## rho_p 0.0622 0.0128
## sigma_x 0.4546 0.0041
## sigma_y 0.467 0.0042
## rho_o -0.5354 1.7234
## tau_x 0.0266 0.0176
## tau_y 0.0288 0.0247
##
## --------------
## YL25
## --------------
## Parameter Estimate Std.Err
## rho_p 0.219 0.0132
## sigma_x 0.6061 0.0054
## sigma_y 0.4665 0.0049
## rho_o -0.9189 0.0762
## tau_x 0.159 0.0317
## tau_y 0.4218 0.0297
##
## --------------
## YL29
## --------------
## Parameter Estimate Std.Err
## rho_p 0.1668 0.0133
## sigma_x 0.5614 0.0054
## sigma_y 0.5372 0.0052
## rho_o -0.8604 0.7478
## tau_x 0.0531 0.0216
## tau_y 0.0492 0.0266
##
## --------------
## YL73
## --------------
## Parameter Estimate Std.Err
## rho_p -0.0781 0.0101
## sigma_x 0.6058 0.0044
## sigma_y 0.5927 0.0043
## . . .
## tau_x 0.0519 0.0178
## tau_y 0.0542 0.0226
##
## --------------
## YL77
## --------------
## Parameter Estimate Std.Err
## . . .
## sigma_x 0.5165 0.0059
## sigma_y 0.464 0.0051
## . . .
## . . .
## . . .
##
## --------------
## YL78
## --------------
## Parameter Estimate Std.Err
## rho_p -0.0461 0.0161
## sigma_x 0.513 0.0058
## sigma_y 0.5159 0.0059
## rho_o -0.6763 0.7341
## tau_x 0.0571 0.0243
## tau_y 0.091 0.0404
extractPredictions <- function(ssm_fit){
predictions <- data.frame(ssm_fit$predicted) |>
sf::st_as_sf() |>
sf::st_transform(4326)
coords <- sf::st_coordinates(predictions)
predictions$lon <- coords[, 1]
predictions$lat <- coords[, 2]
new_data <- predictions |>
data.frame() |>
dplyr::select(id, date, lon, lat)
return(new_data)
}
elk_ssm_list <- elk_ssm_fit$ssm
elk_regdata_list <- lapply(elk_ssm_list,
extractPredictions)
library(ggplot2)
elk_id <- split(elk_mig2, elk_mig2$id)
track_plots <- c()
for(i in c(1:length(elk_id))){
track_plots[[i]] <- ggplot()+
geom_path(data = elk_regdata_list[[i]], aes(x = lon, y = lat), size = 0.7, color = "red")+
geom_point(data = elk_id[[i]], aes(x = lon, y = lat), size=1, color="black")+
theme_classic()+
xlab("Longitude")+ylab("Latitude")+
facet_wrap(~id,scales="free")
}
track_plots
## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
elk_regdata <- do.call("rbind", elk_regdata_list)
library(sf)
elk_reg_utm <- elk_regdata |>
st_as_sf(coords = c("lon","lat"), crs = 4326) |>
st_transform(32611)
coords <- st_coordinates(elk_reg_utm)
elk_regdata$X <- coords[,1]
elk_regdata$Y <- coords[,2]
elk_regdata$Z <- elk_regdata$X + 1i*elk_regdata$Y
getMovementMetrics <- function(dataframe){
move_step <- diff(dataframe$Z)
time_step <- as.numeric(difftime(dataframe$date[-1], dataframe$date[-length(dataframe$date)], "days"))
absolute_turnangle <- Arg(move_step)
relative_turnangle <- diff(absolute_turnangle)
step_length_km <- Mod(move_step)/1000
dataframe$time_step_days <- c(NA, time_step)
dataframe$step_length_km <- c(NA, step_length_km)
dataframe$relative_turnangle <- c(NA, NA, relative_turnangle)
return(dataframe)
}
elk_reg_id <- split(elk_regdata, elk_regdata$id)
elk_reg_id2 <- lapply(elk_reg_id,
getMovementMetrics)
The Lavielle’s method identifies breaking points in a timeseries by using a moving window through the timeseries and a penalized contrast function which determines the optimal number of segments in the timeseries by minimizing the penalized contrast function (Lavielle 2005)
The method examine the variation in the mean and/or variance of a timeseries and identify the optimal number of break points, by fitting the mean, variance or both to the timeseries depending on the number of segments (K). As K increases, the fit increases. However, there should be a number of segments K after which adding more segments does not contributes significantly to improving the fit. The penalized contrast function adjust the trade-off between fit and K
library(adehabitatLT)
str(elk_reg_id2[["YL73"]])
## 'data.frame': 3076 obs. of 10 variables:
## $ id : chr "YL73" "YL73" "YL73" "YL73" ...
## $ date : POSIXct, format: "2004-02-20 00:00:00" "2004-02-20 02:00:00" ...
## $ lon : num -116 -116 -116 -116 -116 ...
## $ lat : num 51.8 51.7 51.7 51.7 51.7 ...
## $ X : num 602201 599949 599230 599220 599239 ...
## $ Y : num 5734480 5733904 5733606 5733601 5733364 ...
## $ Z : cplx 6e+05+5.73e+06i 6e+05+5.73e+06i 6e+05+5.73e+06i ...
## $ time_step_days : num NA 2 2 2 2 2 2 2 2 2 ...
## $ step_length_km : num NA 2.3245 0.778 0.0113 0.2374 ...
## $ relative_turnangle: num NA NA 0.1424 0.0303 1.2283 ...
ggplot(data = elk_reg_id2[["YL73"]], aes(x = date, y = lat))+
geom_path(color="darkgrey") +
geom_point(size=0.5) +
theme_classic() +
xlab("Date") + ylab("Latitude")
decide a minimum number of relocations within a segment (Lmin), the maximum number of segments we want the function to try (Kmax) and if it should look at the difference in mean, variance, or both (mean, var, meanvar)
elk_example_data <- elk_reg_id2[["YL73"]]
elk_lavielle_example <- lavielle(elk_example_data$lat, Lmin = 50, Kmax = 10, type = "mean")
chooseseg(elk_lavielle_example)
## K D Jk
## 1 1 Inf 3075.00000
## 2 2 4.914028021 860.63701
## 3 3 1.313698068 282.70349
## 4 4 0.379474355 142.24697
## 5 5 0.000609291 128.15990
## 6 6 0.001766675 114.27573
## 7 7 0.021777072 100.97988
## 8 8 -0.014404272 94.93606
## 9 9 0.013929530 84.09544
elk_lavielle_example_path <- findpath(elk_lavielle_example, 3)
breakpoints <- unlist(elk_lavielle_example_path)
breakpoints
## [1] 1 892 893 2924 2925 3076
breakpoints_df <- data.frame(Row_ID = breakpoints, segment = c(1,1,2,2,3,3))
elk_example_data$Row_ID <- cumsum(rep(1, nrow(elk_example_data)))
elk_example_data_lav <- merge(elk_example_data, breakpoints_df, by="Row_ID", all.x=TRUE)
elk_example_data_lav <- tidyr::fill(elk_example_data_lav, segment, .direction = 'down')
head(elk_example_data_lav)
## Row_ID id date lon lat X Y
## 1 1 YL73 2004-02-20 00:00:00 -115.5194 51.75189 602200.8 5734480
## 2 2 YL73 2004-02-20 02:00:00 -115.5522 51.74711 599949.0 5733904
## 3 3 YL73 2004-02-20 04:00:00 -115.5627 51.74456 599230.3 5733606
## 4 4 YL73 2004-02-20 06:00:00 -115.5628 51.74452 599220.1 5733601
## 5 5 YL73 2004-02-20 08:00:00 -115.5626 51.74239 599239.2 5733364
## 6 6 YL73 2004-02-20 10:00:00 -115.5575 51.73933 599602.2 5733031
## Z time_step_days step_length_km relative_turnangle segment
## 1 602201+5734480i NA NA NA 1
## 2 599949+5733904i 2 2.32454229 NA 1
## 3 599230+5733606i 2 0.77796897 0.14237354 1
## 4 599220+5733601i 2 0.01127217 0.03027534 1
## 5 599239+5733364i 2 0.23737279 1.22829486 1
## 6 599602+5733031i 2 0.49274756 0.74708438 1
The first passage time (FPT) is a parameter often used to describe the scale at which patterns occur in a trajectory. For a given scale r, it is defined as the time required by the animals to pass through a circle of radius r.
Fauchald & Tveraa (2003) proposed another use of the FPT. Instead of computing the mean of FPT, they propose the use of the variance of the log(FPT). This variance should be high for scales at which patterns occur in the trajectory (e.g. area restricted search). This method is often used to determine the scale at which an animal seaches for food
library(mapview)
elk_example_track <- elk_example_data |>
st_as_sf(coords=c("lon","lat"), crs= 4326) |>
st_union() |>
st_cast("LINESTRING")
mapview(elk_example_track)
transform the data into ltraj object
elk_example_sp <- elk_example_data |>
st_as_sf(coords=c("lon","lat"), crs= 4326) |>
st_transform(32611) |>
st_cast("POINT") %>% as_Spatial()
elk_example_traj <- as.ltraj(coordinates(elk_example_sp), date=elk_example_data$date, id = elk_example_data$id)
summary(elk_example_traj)
## id burst nb.reloc NAs date.begin date.end
## 1 YL73 YL73 3076 0 2004-02-20 2004-11-02 06:00:00
head(elk_example_traj[[1]])
## x y date dx dy
## YL73.X 602200.8 5734480 2004-02-20 00:00:00 -2251.83611 -576.828406
## YL73.X.3 599949.0 5733904 2004-02-20 02:00:00 -718.61809 -298.033141
## YL73.X.6 599230.3 5733606 2004-02-20 04:00:00 -10.27673 -4.631476
## YL73.X.9 599220.1 5733601 2004-02-20 06:00:00 19.18594 -236.596153
## YL73.X.12 599239.2 5733364 2004-02-20 08:00:00 362.94768 -333.270359
## YL73.X.15 599602.2 5733031 2004-02-20 10:00:00 116.54171 8.977010
## dist dt R2n abs.angle rel.angle
## YL73.X 2324.54229 7200 0 -2.8908256 NA
## YL73.X.3 777.96897 7200 5403497 -2.7484521 0.14237354
## YL73.X.6 11.27217 7200 9588981 -2.7181767 0.03027534
## YL73.X.9 237.37279 7200 9658265 -1.4898819 1.22829486
## YL73.X.12 492.74756 7200 10016404 -0.7427975 0.74708438
## YL73.X.15 116.88694 7200 8853351 0.0768765 0.81967397
plot(elk_example_traj)
radii - numeric vector for the radii of the circles
units - time units of the results
elk_example_fpt <- fpt(elk_example_traj, radii = seq(50, 10000), units="hours")
peak in variance at spatial scale of interest
varlogfpt(elk_example_fpt, graph=TRUE)
elk_example_fpt <- fpt(elk_example_traj, radii = 10000, units="hours")
plot(elk_example_fpt, scale = 9000, warn = FALSE)
elk_example_data$FPT <- elk_example_fpt[[1]]$r1
ggplot() +
geom_path(data = elk_example_data, aes(x = date, y = scale(FPT)), color="darkgrey", size = 1) +
geom_path(data = elk_example_data, aes(x = date, y = scale(lat)), color ="orange", alpha = 0.6, size = 1) +
geom_point(data = elk_example_data, aes(x = date, y = scale(FPT))) +
theme_classic() +
ylab("Scaled FPT (grey) vs Latitude (orange)")
The “smoove” package allows for multiple continuous-time movement models (correlated velocity models or CVMs) to be fit to irregular tracking data to identify “modes” of behavior. Each model comes with it’s own set of parameters and best models are determined using AIC/BIC and maximum likelihood. CVM model options include: Uncorrelated CVM (UCVM), Advective CVM (AVCM), Rotational CVM (RCVM), and Rotational-Advective CVM (RACVM). A single-change point can also be identified assuming just a UCVM process and sudden changes in the UCVM parameters, time scale, and root mean square speeds.
library(devtools)
install_github("EliGurarie/smoove")
library(smoove)
elk_utm <- elk_regdata |>
st_as_sf(coords = c("lon","lat"), crs = 4326) |>
st_transform(32611)
coords <- st_coordinates(elk_utm)
elk_regdata$X <- coords[,1]
elk_regdata$Y <- coords[,2]
elk_regdata$Z <- elk_regdata$X + 1i*elk_regdata$Y
elk_example <- elk_regdata |>
dplyr::select(id, date, X, Y, Z) |>
dplyr::filter(id=="YL73")
with(elk_example, scan_track(x = X, y = Y, main = "YL73"))
specify window size (how long does behavior last?)
try: 200 days
takes awhile to run
elk_sweep <- with(elk_example, sweepRACVM(Z = Z, T = date, windowsize = 200, time.unit = "days", windowstep = 30, model = "RACVM", progress = TRUE))
parallel example
library(doParallel)
cl <- parallel::makeCluster(parallel::detectCores())
doParallel::registerDoParallel(cl)
elk_sweep <- with(elk_example, sweepRACVM(Z = Z, T = date, windowsize = 200, time.unit = "days", windowstep = 30, model = "RACVM", progress = TRUE, .parallel = TRUE))
colors represent the relative log-likelihood profile for a single window
looking for distinct peaks (significant change points)
plotWindowSweep(elk_sweep, main = "YL73")
Can specify the clusterwidth here- I used widths ranging from 1-6, based on warnings given with fitting the function for each id
I used all models as candidates (UCVM, ACVM, RCVM, RACVM)
I additionally used AIC to determine models, as BIC can be overly conservative
elk_cp <- findCandidateChangePoints(windowsweep = elk_sweep, clusterwidth = 6) |>
getCPtable(modelset = "all",criterion="AIC")
estimate phases
elk_phases <- elk_cp |>
estimatePhases()
## phase start end model eta tau rms
## 1 I 2004-02-20 00:00:00 2004-05-05 10:00:00 UCVM 10747.467 0 10747.467
## 2 II 2004-05-05 10:00:00 2004-11-02 06:00:00 UCVM 8627.125 0 8627.125
elk_phases_summ <- summarizePhases(elk_phases)
str(elk_phases_summ)
## 'data.frame': 2 obs. of 7 variables:
## $ phase: chr "I" "II"
## $ start: POSIXct, format: "2004-02-20 00:00:00" "2004-05-05 10:00:00"
## $ end : POSIXct, format: "2004-05-05 10:00:00" "2004-11-02 06:00:00"
## $ model: chr "UCVM" "UCVM"
## $ eta : num 10747 8627
## $ tau : num 0 0
## $ rms : num 10747 8627
elk_phases_data <- merge(elk_example, elk_phases_summ, all = TRUE, by.x = c("date"), by.y = c("end")) |> arrange(date) |>
tidyr::fill(phase, start, model, eta, tau, rms, .direction = "updown")
str(elk_phases_data)
## 'data.frame': 3077 obs. of 11 variables:
## $ date : POSIXct, format: "2004-02-20 00:00:00" "2004-02-20 02:00:00" ...
## $ id : chr "YL73" "YL73" "YL73" "YL73" ...
## $ X : num 602201 599949 599230 599220 599239 ...
## $ Y : num 5734480 5733904 5733606 5733601 5733364 ...
## $ Z : cplx 6e+05+5.73e+06i 6e+05+5.73e+06i 6e+05+5.73e+06i ...
## $ phase: chr "I" "I" "I" "I" ...
## $ start: POSIXct, format: "2004-02-20 00:00:00" "2004-02-20 00:00:00" ...
## $ model: chr "UCVM" "UCVM" "UCVM" "UCVM" ...
## $ eta : num 10747 10747 10747 10747 10747 ...
## $ tau : num 0 0 0 0 0 0 0 0 0 0 ...
## $ rms : num 10747 10747 10747 10747 10747 ...
ggplot(data = elk_phases_data, aes(x = date, y = Y, color = model))+
geom_path(size=1.5)+
xlab("Datetime")+ylab("Latitude")+
theme_classic()
ggplot(data = elk_phases_data, aes(x = date, y = Y, color = phase))+
geom_path(size=1.5)+
xlab("Datetime")+ylab("Latitude")+
theme_classic()